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We study the rotation curves of ultralight BEC dark matter halos. These halos are long lived 
solutions of initially rotating BEC fluctuations. In order to study the implications of the rotation 
characterizing these long-lived configurations we consider the particular case of a boson mass m = 
lCU 23 eV/c 2 and no self-interaction. We find that these halos successfully fit samples of rotation 
curves (RCs) of LSB galaxies. 


PACS numbers: 95.35.+d,98.62.Gq 

I. INTRODUCTION 

The idea of an ultralight spinless particle as dark 
matter originates at cosmic scales, when it was discov¬ 
ered that the predicted mass power spectrum of fluctua¬ 
tions made of a free real scalar field, minimally coupled 
to General Relativity in an affective theory, with mass 
TO^ ~ 1CU 22 — 10 _23 eU could ameliorate the problem 
of small structures associated to the standard cold dark 
matter model [Tf . An interesting property of interpreting 
this scalar field as a spinless boson is that its condensa¬ 
tion temperature is T c ~ l/m 5 / 3 ~ TeV for ~ 10 -22 . 
This interpretation permits the formulation of the struc¬ 
ture formation problem as an evolution problem ruled 
by the evolution of the condensate, namely the Gross- 
Pitaevskii equation Q that describes the evolution of 
the condensate 0: with a potential due to the gravity 
of the fluctuation itself. It was until very recently that 
the dynamic structure formation problem was analyzed 
for such a condensate for the case of an ultralight boson 
with no self-interaction Q. 

Aside of the different emerging questions on the anal¬ 
ysis of structure formation within the model, the prob¬ 
lem at local scale remains. This problem involves the 
analysis of the collapse, formation, evolution and virial- 
ization of BEC halos that would be ruled by the evo¬ 
lution of a BEC trapped by its own gravitational field, 
that is, the Gross-Pitaevskii equation coupled to Pois¬ 
son equation, the GPP system. In fact this local prob¬ 
lem has been studied from different angles and some ad¬ 
vances have been achieved, which are summarized briefly 
as follows. Historically, the first time that a solution to 
the GPP system in spherical symmetry was presented, it 
was in one of the original papers defining Boson Stars, 
because it happens that the low-energy and Newtonian 
limit of the Einstein-Klein-Gordon for a complex scalar 
field is exactly the GPP system of equations, and nu¬ 
merical solutions of stationary equilibrium configurations 
were constructed Q. Later on, other equilibrium config¬ 
urations -dubbed gravitational atoms- were constructed, 
which were solutions to the same GPP system, however 
for wave functions with nodes; the fact that the solu¬ 


tion wave functions had nodes invited to call them ex¬ 
cited state solutions as an analogy to the solution of the 
Schrodinger equation for the Hydrogen atom [gi, 0|. The 
excited state solutions were in fact applied as models of 
galactic dark matter, because the RCs of the resulting 
configurations had very similar properties of observed 
galactic RCs @. These excited states were shown to be 
unstable though @,[1]. At the same time, other flavors or 
approaches were developed, for instance another version 
of the model using the GPP system to explain the RCs 
was dubbed the model of quintessential halos 00; 
another example is the fuzzy dark matter model, that is 
actually a ID version (not even spherically symmetric in 
3D) of the analysis of the collapse of fluctuations obeying 
the GPP system 0 - A very complete analysis of BEC 
spherical solutions can be found in i ;I Li 

Some other BEC dark matter halo models include 
those constructed in the Thomas-Fermi limit (TF), in 
which the self-interaction among bosons dominates over 
the kinetic term in Schrodinger equation 0. Finite tem¬ 
perature effects have also been considered in this regime 
0, and it has been shown to be a possible solution of 
the galactic cusp-core problem 0. The life-time of these 
halos has been estimated 0, their energy contents flp 
and the condensation process has also been studied 0]. 
Recently, the full collapse process of BEC dark matter 
fluctuations obeying the GPP system has been analyzed 
in 0 , [22| . When compared with observations, galac¬ 
tic halos constructed in the TF limit are capable of fit¬ 
ting RCs, although some uncertainty has been mentioned 
with respect to the possibility of their formation [0 . 

In this paper we do not consider the model in the 
TF limit. Instead we consider the system is not dom¬ 
inated by the self-interaction nor the kinetic contribu¬ 
tion and we assume the whole GPP system to hold. As 
mentioned above, equilibrium configurations can be con¬ 
structed, however some steps forward have been done 
in more dynamically general situations. For instance, it 
has been shown that initial fluctuations of the density 
of the condensate quickly virialize after the turnaround 
point, in which the cosmic expansion ceases to dominate 
over the self-gravity of the structure 0], through the 
gravitational cooling HU. It was also shown that in 
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general the evolution of an initial spherical fluctuation 
ends up approaching one of the ground state equilibrium 
solutions [9( and also that non-spherical initial fluctua¬ 
tions approach these equilibrium solutions as well [25j . 
Thus equilibrium configurations are stable under a va¬ 
riety of perturbations, including non-spherical perturba¬ 
tions, and also show an attractor type of behavior. These 
properties together are the reason to consider ground 
state solutions a very appealing candidate to be dark 
matter halos. 

In this general scenario the dynamical properties of 
the halo formation and virialization has therefore been 
explored, however, unlike the models in the TF limit, 
RCs have not been explored for already relaxed config¬ 
urations, and the goal of the present paper is precisely 
to show their potential at fitting RCs. At this point two 
alternatives have been explored so far. On the one hand, 
even though excited states are unstable, the superposi¬ 
tion of ground and excited states turned out to be sta¬ 
ble and show appealing galactic RCs [26|. On the other 
hand, some rotation has been added to ground state con¬ 
figurations that disperse away the density of the conden¬ 
sate and show also appealing RCs [27]]. The possibility 
that BEC dark matter halos may have rotation has been 
pointed out before in [2S|, where spheroid and ellipsoid 
analytic solutions to the GPP system with rotation are 
studied as rotating BEC dark matter halos in various 
scenarios, and the results particularly focus on the pos¬ 
sibility of vortex formation, but little is said about the 
comparison with observations, say RCs. 

What we do in this paper is to try to answer the ques¬ 
tion of whether the long-lived rotating or spherical con¬ 
figurations in [23] are capable of fitting observed RCs. 
For this we focus on a sample of dark matter dominated 
LSB galaxies, for which, as a first approximation, lu¬ 
minous matter would not contribute significantly to the 
total mass of the system and thus behave as test parti¬ 
cles. As a workhorse we also consider the BEC without 
self-interaction, in order to be consistent with the only 
structure formation analysis so far with BEC dark matter 
at cosmological scale and use the same ultralight boson 
mass 10 -23 eV/c 2 Q. 

The paper is organized as follows, in section [TT| we 
present the conventions and numerical methods used to 
solve the GPP system and the methods used to calculate 
the RCs. In section eh we fit actual RCs and in IIVI we 
discuss the results. 


II. THE GPP SYSTEM 


A. Equations and numerical scaling 


Units and scaling. The Gross-Pitaevskii-Poisson 
(GPP) system of equations describing the evolution of 
a Bose condensate is 


h' 2 - , - 

lh—=r = - v 2 t 

dt 2m 

V 2 P = 4ttGto|T| 2 , 


■V 'F 


2nh 2 c 


|4d 2 4/, 


(1) 


where in general the wave function depends on space and 
time 4' = 41 (t, x), m is the mass of the boson, V is the 
gravitational potential acting as the condensate trap, a 
is the scattering length of the bosons. This is a coupled 
system of an evolution equation for 41 with a potential 
that is solution of Poisson equation sourced by |4/| 2 . 

Before integrating © it is important to remove the 
constants using the following change of variables 4* = 

V4irGh ,f, ™ _ mc~ - _ me ~ i _ me; 1 _ me 2 X 
me 2 x — h x ’ y ~ h y> z — h z ’ 1 — h b 

V = a —> so the numerical coefficients 

h , h 2 /m, 2i xh 2 /m 2 , AitGm do not appear in ©. Thus 
by fixing the boson mass m the hatted variables are fixed. 


Going further, the system © is invariant under the 
transformation t = A 2 t, x = Xx, y = Xy, z = Xz, 4/ = 
T/A 2 , V = P/A 2 a = A 2 a, for an arbitrary value of 
the parameter A [8J . This rescaling reduces the original 
system © to the following one 


i^- = -iv 2 4> + W + a|4'| 2 4' (2) 

V 2 P = ^l 2 , (3) 


which is the one we solve numerically. The consequences 
of this invariance are extremely relevant for our analysis 
here. The reason is that if one calculates one solution for 
the non-hatted variables, say for A = 1, other solutions 
are automatically found for any other value of A, that is, 
a new solution for the hatted and tilde variables is found 
just by using a new value of A. 


Evolution. The solution is calculated considering 
Schrodinger equation as an evolution equation for 4 1 , 
and Poisson equation as a constraint that has to be 
solved every time it is required during the integration 
of Schrodinger equation. 

We approximate the GPP system m using finite 
differences on a uniformly discretized grid on a spa¬ 
tial domain described with cartesian coordinates. We 
solve the system on the spatial domain [x m i n , x max ] x 
[ymin^ymax] x \zminiZmax\ using the 3D Fixed Mesh Re¬ 
finement (FMR) code developed by us and tested in 127 1. 
The evolution uses a method of lines with a Runge-Kutta 
integrator. On the other hand, Poisson equation is solved 
reducing the problem to a 2D slice on the plane x + y = 0 
because we only deal with axial configurations in this 
work; then the solution is interpolated back into the 3D 
mesh. The algorithm used to solve Poisson equation on 
the 2D slice is a successive overelaxation method with 
optimal acceleration parameter. 
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B. Rotation curves 

In order to construct the rotation curve of one of our 
BEC halos, we place various detectors at which we cal¬ 
culate the tangential velocity v of a test particle. The 
detectors are located at a set of points along the x- 
axis. We assume the test particle describes a circular 
orbit which implies that the tan gent velocity of the test 
particle is v(r) = y/2GM(r)/r, where r is the distance 
from the coordinate center to the detector, and M(r) 
is the mass contained within a sphere of radius r, and 
calculated as a volume integral in the 3D domain as 
M = f \^\ 2 dxdydz. Explicitly, if a particular detector 
is located at (x, y, z) = (x d , 0,0) we calculate v(xd) as 

v 2 {x d ) = [ |^| 2 dxdydz, (4) 

\Xd\ J 

where the volume integral is calculated within a sphere 
of radius x d on our 3D cartesian grid. 


C. Initial data 


Our BEC halos are the relaxed long-lived configura¬ 
tions resulting from the evolution of an initial fluctuation 
we call protohalo. 

The initial data we choose for a proto-halo, is a ground 
state equilibrium configuration added with an initial an¬ 
gular momentum. When the added angular momentum 
is zero, the protohalo is actually a ground state equilib¬ 
rium configuration stable and long-lived since the begin- 
ing. When some angular momentum is added initially 
to an equilibrium configuration, it is not a ground state 
equilibrium configuration anymore, and it takes a time 
for the configuration to relax and become a long-lived 
structure that plays the role of one of our BEC halos. 
Thus the construction of initial data requires a descrip¬ 
tion of the equilibrium configurations. 

Equilibrium configurations. These are built assuming 
the wave function is spherically symmetric, a reason to 
use spherical coordinates, and it is also assumed it is har¬ 
monically time dependent 4/ = vf (r, t) = e lut ip(r), which 
immediately guarantees that |'k| 2 is time-independent 
and therefore the potential V = V(r) as well. Under 
these conditions the GPP system ([3]) reduces to 


d 2 ip 2 dip 


dr 2 
d 2 V 

dr 2 


r dr 
2 dV 
r dr 


2 (V + a\ip\ 2 + uj)ip, 
4tt|0| 2 , 


(5) 


which is an eigenvalue problem for 0(r) with the bound¬ 
ary conditions of wave function smoothness at the ori¬ 
gin and isolation at infinity, that is 0(r —> oo) —> 0. 
This problem is then solved numerically in a finite spa¬ 
tial domain. The boundary condition for the gravi¬ 
tational potential is V = -M/r for r oo where 


M = f \ip(r)\ 2 r 2 dr as described in Q. The solution is 
such that given a value for "0(0) there is a unique eigen¬ 
value ui satisfying the boundary conditions and plays the 
role of the eigen-energy of the system. 

The solutions of the eigenvalue problem © can have 
nodes or not. Those configurations with nodes are called 
excited state solutions whereas those without nodes are 
called ground state solutions. In Q was shown in detail, 
using numerical simulations, that excited state configu¬ 
rations are unstable and in fact decay into one of the 
equilibrium configurations through a relaxation process 
called gravitational cooling that consists in the ejection 
of density of probability to infinity am. 

With this in mind we are only interested in ground 
state configurations. We then remark that there is a 
configuration for each different value of 0(0) that in the 
end has a specific associated value of ui. Two comments 
are in turn: i) given the rescaling property of the GPP 
system, it is not necessary to construct all the possible 
wave functions and density profiles for each value of 0(0), 
because given one, say 0(0) = 1, one can construct all the 
other possible ground state configurations using simply 
different values of the scale parameter A; ii) since we are 
interested in the density of the BEC, that is |T| 2 , uj does 
not play a relevant role in our analysis and in fact -if 
required- the value of this frequency is associated to a 
specific central value of the wave function 0(0). 

As example solution of (5j), we show in Fig. [Qthe den¬ 
sity profile and the rotation curve corresponding to the 
case 0(0) = 1 for a = 0 and a = 1. The introduction of 
a positive self-interaction allows more massive and wider 
BEC dark matter distributions as shown in the Figure. 
One of the reasons to introduce a self-interaction term, is 
that it might imply better RCs, however, as can also be 
seen in Fig. |T] the rotation curves with or without self¬ 
interaction show a pretty similar shape. Furthermore, 
the fact that self-interaction produces wider configura¬ 
tions does not mean that the density of the BEC is more 
disperse in space, instead it is more compact. Consid¬ 
ering the examples in the figure, we have that for the 
case of a = 0 the mass and compactness are M = 25.91, 
M /rgs = 6.61, whereas for a = 1, these quantities are 
M = 60.79 and M/rgs = 13.94. Here rg 5 is the ra¬ 
dius containing 95% of the total mass of the configura¬ 
tion. The failure at improving RCs when adding self¬ 
interaction, is one of the motivations to introduce some 
rotation, because rotation actually disperses away some 
of the density concentrated in the center, and the reason 
why we set a = 0 from now on. 

Addition of rotation. Our code is three dimensional, 
and the eigenvalue problem we solve to construct spheri¬ 
cally symmetric configurations (JS]| uses spherical coordi¬ 
nates. Thus, in order to set initial data we interpolate the 
solution 0(r) of (JSJ) into our 3D numerical grid and ob¬ 
tain 4>(a;, y,z,t = 0). Once this has been done, we apply 
a rotation by redefining the wave function 4/ = e _lL ' n0 'I'. 
In our particular case of rotation around the z-axis we 
use 9 = arctan {y/x) with L = L z z. We parametrize the 
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FIG. 1: Properties of the ground state spherically symmetric 
configuration for the central wave function value ^i(O) = 1. In 
the top panel we show the density of the BEC and the effect 
of the self-interaction term. In the bottom panel we show the 
effects on the rotation curve with and without self-interaction. 


rotation by choosing L z = xp y — yp x to be a constant. 


D. Simulations in code units 

The criterion to consider a configuration as long-lived 
halo is that the rotation curve remains nearly time- 
independent after a transient time. We show in Fig. [2] 
the RC evolution of four configurations for A = 1 and 
L z = 0 corresponding to formal equilibrium and three 
other values L z = 0.82,0.85, 0.87. In the case L z = 0 the 
system remains time independent as expected, because it 
is a spherically symmetric equilibrium configuration; in 
the other three cases the RCs start high and while the 
density disperses away the RC flattens and starts becom¬ 
ing time independent. It can be seen how the snapshots 
start packing around a nearly time-independent profile. 
Such long-lived profile is the one we choose to fit the 
observed galactic RCs in this paper. The behavior for 
other values 0.82 < L z < 0.87 should show a similar be¬ 
havior bounded by these two values, however the use of 
the three chosen values suffice to illustrate the capability 
of the rotating BEC halo model to fit RCs. 

Additionally, we made sure that the total energy of 
the system is negative and the gravitational potential 
depth approaches an asymptotic long-living behavior in 




FIG. 2: Snapshots at various times of the RCs for four val¬ 
ues of L z , which stabilize around a given profile. The time 
direction of the snapshots is from top to bottom in all the 
cases. 

the scale of Gyr. 

E. From code units to physical units 

The code and physical units are related once the scale 
invariance parameter A is fixed, which we do by choosing 
the spatial coordinate units. Specifically, since A = 

(or equivalently the spherical coordinate r), if physi¬ 
cal measurements of space x,y,z are given in kpc, we 
choose the code coordinates x, y , £ to represent kpc as 
well. Thus it suffices to write the factor — in kpc. For 
m = 10 - 23 eV/c 2 its value is A := A 0 = ^[ k P c ]j^ = 
0.0006389. 

However, the results obtained with the code units are 
valid for any value of A, and in fact we parametrize the 
value of A = aAo with a a regulation parameter that 
shifts the relation between code and physical units in 
such a way that if a > 1 a code unit of length represents 
more than a kpc, whereas a < 1 does the opposite. 

Once the value of A is fixed, we show how the rele¬ 
vant quantities translate from the code units into physi¬ 
cal units. 

Rotation Curve velocity. We start by writing in 
physical units and substitute the scaling transformations 
for tilde and hat variables 

-2 2 G /■,? |2 A 2 c 2 2 G 

v = 7 — / OF \ dxdydz = — — 7—7 / OF \ dxdydz, 

\Xd | J 27T \x d \ J 

and then, writing c in units of km/s and A = Aon = 
0.0006389a, the velocity v in km/s in terms of the veloc¬ 
ity in code units v relate through 

, c[km/sl , cfkm/sl 

v = A— _ v = aAo — _ v = 76.597au[km/s]. 

v27t v27t 
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FIG. 3: Variety of RCs that can be generated for different 
values of L z . Each curve of each plot corresponds to a differ¬ 
ent value of a. We show the curves generated for a > 1 only, 
and the bold line corresponds to a = 1. 


Length. Even though we already fixed the length units 
we want to specify the transformation x = j = 
Coordinates y, z scale identically. 

The scaling of velocity and length tells us about the 
bundle of RCs that can be constructed. By choosing 
a > 1 the RC in physical units shows a higher rotation 
curve, but at the same time it shortens the size of a halo 
and conversely for 0 < a < 1. In Fig. [3] we show the 
bundle of different RCs produced by different values of 
a and the four chosen values of L~. The larger the value 
of L z the more disperse the configuration is and smaller 
the velocity of test particles. 

With these specifications, the simulations we used to 
construct the long-lived configurations where carried out 
in a domain of size [—80kpc, 80kpc] 3 for a = 1 and 
four refinement levels with resolutions 1.25kpc cover¬ 
ing the domain [—80kpc, 80kpc] 3 , 0.625kpc covering the 
domain [—40kpc, 40kpc] 3 , 0.3125kpc covering the do¬ 
main [—20kpc, 20kpc] 3 , 0.15625kpc covering the domain 
[—lOkpc, lOkpc] 3 . This domain was sufficient for the con¬ 
figurations to relax and approach a long-lived state. 


Galaxy 

Best-Lz 

a 

ES03020120 

0.82 

1.7 

ES03050090 

0 

0.325 

ES04880049 

0 

0.45 

U4115 

0 

0.53 

U11557 

0.87 

3.2 

U11616 

0.82 

3.5 


TABLE I: Best values of L z and a that fit the RCs. 


by 0 


vrm{t) 



where R indicates the radius at which the density of dark 
matter is zero. We have used the parameters R : p 0 found 
in 0 to produce their plots. Also for comparison and 
control we have included the profile due to a Pseudo 
Isothermal profile (PI) which rotation curve profile is 
given by [3l| and later used to directly fit observations 


(W^larcan (j_)). (7) 

We use the best fitting parameters p 0 ,Rpj found in [30| . 
We have extended the spatial domain on purpose, so that 
future experimental points may decide between the mod¬ 
els constructed assuming the Thomas-Fermi limit and 
our more general model. 

A clear difference between our halo model of BEC dark 
matter and that assuming the Thomas-Fermi limit, is 
that in the later case there is a cut off of the density 
at a finite radius that makes the rotation curve to drop 
drastically. Clear cases are the galaxies ES03020120, 
ES04880049 and U11616 in Fig. [J1 where the rotation 
curve of the RM model fits only in the inner parts, but 
clearly drops further out. In these cases our rotating 
model performs much better and shows a trend more 
similar to the PI control profile. 


vpi(r) = i AnGpoRjz 


pi 


III. FITTING ROTATION CURVES 

As mentioned before, we use values of L z = 
0, 0.82, 0.85, 0.87, that we found empirically to allow 
long-lived configurations. For each of these values we 
track a value of A that fits the RCs data. The sample we 
use to fit the RCs of our rotating halos is a subsample of 
LSB galaxies without luminous components in (29| . 

We were able to successfully fit a broad sample of LSB 
galaxies, as shown in Fig. []] and the best fitting pa¬ 
rameters are shown in Tabic [I] For comparison we have 
also plotted the RCs due to halos made of 10 _22 eV/c 2 
bosons in the Thomas-Fermi limit constructed by Robles 
& Matos (RM) in [30[, which have a velocity profile given 


IV. DISCUSSION AND CONCLUSIONS 

Given the diversity of RCs in LSB galaxies, we explored 
the chances that different values of L z could also explain 
such diversity. We have shown that fixing m and a (in 
our case a = 0), a bundle of RCs with different profiles 
can be constructed for each value of L z . 

A more general analysis of the parameter space in¬ 
cludes the variation m and a, which would imply a three 
dimensional parameter space. This situation rather de¬ 
fines an inverse problem. That is, an analysis looking to¬ 
ward the BEC dark matter model will require data analy¬ 
sis of a universal sample of galaxies in order to fix a single 
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FIG. 4: RC for the galaxies ES03020120, ES03050090, ES04880049, U4115, U11557 and U11616. We extend the domain in 
order to show the drastic fall of the curve with halos constructed in the TF limit. 


value of m and a, so that the same bosons are the same 
dark matter in every halo, a mistake usually overseen in 
models in the Thomas-Fermi limit, where different val¬ 
ues of a are found for different galaxies, which is somehow 
inconsistent (e. g. 0 Jill )• 

Our results show that rotating BEC dark matter halos, 
and non-rotating ground state equilibrium configurations 
are an option worth to study in more general samples of 
galaxies, in which luminous matter has a more consid¬ 
erable contribution. This however will require the de¬ 
velopment of a code that in the minimal case solves the 
GPP system coupled to Euler equations to describe the 


luminous matter. 

Finally, even though our working hypotheses of an ul¬ 
tralight boson and zero self-interaction are similar to the 
only structure formation analysis so far in Q, it would be 
also interesting to consider the recent restrictions found, 
imposed by BBN on m and a [32], H3 • 

Acknowledgments 

We appreciate the comments from the anonymous ref¬ 
eree. This research is partly supported by grants: CIC- 




























7 


UMSNH-4.9 and CONACyT 106466. F.D.L-C grate¬ 
fully acknowledges a postdoctoral DGAPA-UNAM grant. 
F.S.G. acknowledges support from the CONACyT pro¬ 


gram for sabbatical visits in foreign countries and the 
kind hospitality of the Physics and Astronomy Depart¬ 
ment at UBC, where part of this paper was developed. 


[1] V. Sahni and L. Wang, Phys. Rev. D 62, 103517 (2000). 
T. Matos and L. A. Urena-Lopez, Phys. Rev. D 63, 
063506 (2001) 

[2] E. P. Gross, Ann of Phys 4 (1958) 57. E. P. Gross, J. 
Math. Phys. 4 (1963) 195. L. P. Pitaevskii, Sov. Phys. 
JETP 13 (1961) 451. 

[3] L. Arturo Urena-Lopez Phys. Rev. D 90, 027306 (2014). 

[4] Hsi-Yu Schive, Tzihong Chiueh and Tom Broadhurst, 
Nature Physics 10, 496499 (2014). 

[5] R. Rufiini and S. Bonazolla, Phys. Rev. 187 (1969) 1767. 

[6] S-J Sin, Phys. Rev. D 50 (1994) 3650. S.U. Ji and S.J. 
Sin, Phys. Rev. D 50 (1994) 3655. 

[7] R. Ferrell and M. Gleiser, Phys. Rev. D 40 (1989) 2524. 

[8] F. S. Guzman and L. A. Urena-Lopez, Phys. Rev. D 69 
(2004) 124033. arXiv:gr-qc/0404014 

[9] F. S. Guzman and L. A. Urena-Lopez, ApJ. 645 (2006) 
814-819. 

[10] A. Arbey, J. Lesgourgues, P. Salati, Phys. Rev. D 
64,123528 (2001). 

[11] A. Arbey, J. Lesgourgues, and P. Salati, Phys. Rev. D 
68, 023511 (2003). 

[12] W. Hu, R. Barkana, and A. Gruzinov, Phys. Rev. Lett. 
85, 1158 (2000). 

[13] P. H. Chavanis, L. Delfini, Phys. Rev. D 84, 043532 

( 2011 ) 

[14] P. H. Chavanis, Phys. Rev. D 84, 043531 (2011) 

[15] C. G. Boehmer and T. Harko, JCAP 07 (2007) 025. 

[16] Tiberiu Harko, Eniko J. M. Madarassy, JCAP 01 (2012) 

020 . 


[17] T. Harko JCAP 1105:022,2011 

[18] F. S. Guzman, F. D. Lora-Clavijo, J. J. Gonzalez-Aviles, 
F. J. Rivera-Paleo, JCAP 12 (2013) 015. 

[19] J. C. C. de Souza, M. O. C. Pires JCAP 03 (2014) 010 

[20] T. Harko, Phys.Rev.D83:123515,2011. 

[21] T. Harko, Phys. Rev. D 89, 084040, 2014 

[22] P. H. Chavanis, A & A 537, A127 (2012) 

[23] F. S. Guzman and L. A. Urena-Lopez, Phys. Rev. D 68 
(2003) 024023.. 

[24] E. Seidel and W-M Suen, Phys. Rev. D 42 (1990) 384. 

[25] A. Bernal and F. S. Guzman, Phys. Rev. D 74 (2006) 
063504. 

[26] L. A. Urena-Lpez, A. Bernal, Phys. Rev. D 82 (2010) 
123535. 

[27] F. S. Guzman, F. D. Lora-Clavijo, J. J. Gonzalez-Aviles, 
F. .J. Rivera-Paleo, Phys. Rev. D. 89 (2014) 063507. 

[28] T. Rindler-Daller and P. R. Shapiro, MNRAS 422 (2012) 
135-161. Ibid. arXiv: 1209.1835. 

[29] W. J. G. de Blok, Stacy S. McGaugh and Vera C. Rubin, 
AJ, 122 (2001) 2396-2427. 

[30] V. H. Robles, T. Matos, MNRAS, 2012, 422, 282-289. 

[31] K.G. Begeman, A. H. Broeils, R. H. Sanders, MNRAS 
249 (1991) 523. 

[32] T. Rindler-Daller, P. R. Shapiro, Mod. Phys. Lett. A 29 
(2014) 1430002. 

[33] Bohua Li, Tanja Rindler-Daller, Paul R. Shapiro, Phys. 
Rev. D 89, 083536 (2014) 



